Green's Function Monte Carlo study of correlation functions in the (2+l)D U(l) 

lattice gauge theory 

C. J. HamerO, R. J. Bursil^ and M. bamaras 
School of Physics, University of New South Wales, Sydney, 2006, Australia. 
^Department of Physics, UMIST, PO Box 88, Manchester, M60 1QD, UK. 

A "forward walking" Quantum Monte Carlo (QMC) algorithm has 
been developed to calculate correlation functions for the Hamiltonian 
lattice formulation of U(l) Yang-Mills theory in (2+1) dimensions. It is 
shown that Wilson loops can be calculated with high accuracy. Creutz 
ratios are used to determine the string tension, which agrees with re- 
sults from other approaches. Timelike correlations are used to estimate 
the mass gaps, which agree with series expansion results in the strong 

■ coupling regime. 
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' I. INTRODUCTION 

' The two major variants of lattice gauge theory (LGT) are the "Euclidean" formulation 
Ph ■ of Wilson Q , and the "Hamiltonian" version of Kogut and Susskind Q . In the Euclidean 
regime, classical Monte Carlo simulations have proved to be extremely powerful in extracting 
quantitative predictions from the theory, as first shown by Creutz [||. This approach is 
preferred by an overwhelming majority of lattice gauge theorists at the present time. 

The Hamiltonian formulation is still worthy of study, however. It can provide a valuable 
check of universality, for instance. Lattice gauge theory relies on the fundamental assump- 

■ tion that quantities such as mass ratios calculated in the continuum limit (a critical point of 
fN| ' the lattice model) must be 'universal', i.e. independent of the microscopic lattice structure 

■ or space-time formulation. There is not much real doubt that this is correct, but it is impor- 
tant to provide checks where possible Q . Another reason is that many techniques imported 
from quantum many-body theory and condensed matter physics can be employed in this 
arena, which may give useful results. Examples include the strong-coupling series approach 

'"" j" 1 • ||, the i-expansion ||, the coupled-cluster method and others. Nevertheless, it seems 
O ,' likely that Monte Carlo simulations will provide the most robust and accurate numerical 
D . techniques in this area also. Our aim in this paper is to discuss some further applications 
'"^ ' of these Quantum Monte Carlo methods || . 
^ . The use of quantum Monte Carlo methods in Hamiltonian LGT has a long and somewhat 
chequered history, and lags a good ten years behind the Euclidean developments. The first 
calculations used a strong-coupling basis involving discrete "electric field" link variables, and 
a "Projector Monte Carlo" approach |lO[|llf| , which used the Hamiltonian itself to project 
out the ground state. A later version of this was the "stochastic truncation" approach of 
Allton et al. |[2| . Using this approach one can successfully compute string tensions and mass 
gaps for Abelian models jn|. For non-Abelian models, however, some problems arose [ H . 
Using an electric field representation for the link variables and a Robson- Webber recoupling 
scheme [[l4| at the vertices requires the use of Clebsch-Gordan coefficients or 6j-symbols, 
which are not known to high order for SU(3); and furthermore, the 'minus sign' problem 
rears its head, in that destructive interference occurs between different paths to the same 
final state. It may well be that a better choice of strong-coupling basis, such as the 'loop 
representation', might avoid these problems; but this has not yet been demonstrated. 

Heys and Stump [jl5| and Chin et al. pioneered the use of "Greens Function Monte 
Carlo" (GFMC) or "Diffusion Monte Carlo" techniques in Hamiltonian LGT, in conjunc- 
tion with a weak-coupling representation involving continuous gauge field link variables. 
This was successfully adapted to non-Abelian Yang-Mills theories [ p~T| , with no minus sign 
problem arising. In this representation, however, one is simulating the wave function in 
gauge field configuration space by a discrete ensemble or density of random walkers: it is 
not possible to determine the derivatives of the gauge fields for each configuration, or to 
enforce Gauss's law explicitly, and the ensemble always relaxes back to the ground state 
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sector. Hence one cannot compute the string tensions and mass gaps directly as Hamil- 
tonian eigenvalues corresponding to ground states in different sectors, as one does in the 
strong-coupling representation. Chin, Long and Robson Q thus resorted to "variational 
Monte Carlo" (VMC) techniques to compute mass gaps. They obtained some reasonable 
results; but this approach always suffers from the major drawback that there is an unknown 
systematic error in the results, due to their dependence on the form of the variational wave 
function which is chosen. 

It appears, therefore, that to make unbiased measurements of mass gaps in the weak- 
coupling representation one is forced back to the more laborious approach used in Euclidean 
calculations: namely, to measure an appropriate correlation function, and estimate the mass 
gap as the inverse of the correlation length. In ref. 1T3] the GFMC method was tested on 
the (2+l)D U(l) model using a "secondary amplitude" technique to compute expectation 
values: but this proves to be expensive and prone to bias jl9 20 1. In this paper we will 
show how the standard 'forward-walking' technique used in many-body theory |pl|] can be 
used for this purpose. The forward-walking method has already been applied to lattice spin 
models by Runge |2^] and Samaras and Hamer pofl . Here we apply it to the compact U(l) 
Yang-Mills theory in (2+1) dimensions, which has been a standard test-bed for Hamiltonian 
LGT. 

A brief discussion of the U(l) model is given in Section II. Our Monte Carlo methods are 
outlined in Section III. The GFMC method is briefly summarized, and then the forward- 
walking method for estimating expectation values is discussed, together with a technique for 
measuring timelike correlations. In Section IV the results are presented. Our conclusions 
are summarized in Section V. 



II. MODEL HAMILTONIAN 



In a weak coupling basis the Hamiltonian for the compact U(l) LGT in (2+l)D is given 
by [i§: 

Z ; P 

where Ai is the gauge field variable on link I and 

6 P = A h + A h - A l3 - A u (2) 

is the plaquette variable for a lattice plaquette P, formed by the four links ■ as 
illustrated in Fig. [j](a). We consider a periodic, square lattice of linear dimension L and 
lattice spacing a. The 'strong coupling' parameter x = l/e 4 a 2 = l/g 4 approaches infinity 
in the continuum limit a — > 0. 

This is an interesting model, which possesses some important similarities with QCD (for 
a more extensive review, see for example ref. [p"3[). If one takes the 'naive' continuum 
limit at a fixed energy scale, one regains the simple continuum theory of non-interacting 
photons p5[ ; but if one renormalizes or rescales in the standard way so as to maintain the 
mass gap constant, then one obtains a confining theory of free massive bosons, as discussed 
by Polyakov |2q| , and proven by Gopfert and Mack ^7j. The Hamiltonian version of the 
model has been well studied by a variety of methods: some of the more recent include 
series expansions |28| , f inite-lattice techniques [j30| , the t-expansion pl| , p2f , and coupled- 
cluster techniques |33|-|35|, as well as QMC p^j36| , |37| , |l3| . Quite accurate estimates have 
been obtained for the string tension and mass gaps, which can be used as comparisons for 
our Monte Carlo results. The finite-size scaling properties of the model can be predicted 
using an effective Lagrangian approach combined with a weak-coupling expansion ]2~i| ] , and 
the predictions agree very well with finite-lattice data (l| . 
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III. MONTE CARLO METHODS 



A. Greens Function Monte Carlo 

We use the Green's Function Monte Carlo [GFMC] method |J , which was adapted to the 
U(l) model by Heys and Stump fl5fl , and Chin et al. |l(J. A brief summary of the method 
can be given as follows. 

In a weak-coupling representation, the basis states are taken to be eigenstates of the 
plaquette angles 9p, which can take continuous values. The Hamiltonian (|l]) can be written 
as 

where 

V{0) = -2x E cos °Pi ( 4 ) 
p 

and the plaquette angles Op and link angles Ai are related by equation (Q). The imaginary 
time Schrodinger equation for the system is 

" r) = [- E ^2 + V(Q) E T M&, r), (5) 

i ' 

where Et is a trial energy, representing a constant shift in the zero of energy, which will 
prove useful. The imaginary time evolution operator exp[— (H — Et)t] acts as a projector 
onto the ground state |$ ) : 

|$ > = lim e - T{H - E ^ |* T > (6) 

r — >oo 

for any trial state |^t), provided that |^t) is not orthogonal to |$o)- 

Equation (||) is a diffusion equation in configuration space, and is easily simulated by the 
Green's Function Monte Carlo method. It is assumed that the ground-state wave function 
can be chosen positive everywhere, and it is simulated by the density distribution of an 
ensemble of random walkers {6i} in configuration space, with weights {uii}. The first term 
on the right of Eq. produces diffusion, and is simulated by a Gaussian random walk 
of the members of the ensemble as time proceeds, while the term |V(0) — -Et] produces a 
growth or decay in the density which is simulated by a branching process. 



B. Variational Guidance 



The efficiency and accuracy of the simulation are greatly enhanced by the use of variational 
guidance or importance sampling ||. Let ^x(©) be a variational approximation to the true 
ground-state wave function, and define a new probability distribution 

/(e,r) = $(e,T)* T (e), (7) 

Then the modified imaginary time Schrodinger equation for /(0,r) reads 

- % = - E % + [^( Q ) - E ^f + E ^( (8) 



where 
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is the local energy obtained from the trial function, and 

f (p» 2 ^ T(e) rim 

is a "quantum force" term, which produces a directed drift in the ensemble towards the 
configurations favoured by the trial wave function. By a good choice of "fx and £r the 
"excess local energy" term [_Z?l(6) — £t] can be made very small, which reduces the amount 
of branching necessary, and reduces the statistical fluctuations in the results. 
For small time steps At, an approximate Green's function solution to Eq. (||) is 

G(9 - 6', At) ~ exp{-[£ L (e) - E t ]At} (11) 
/ 1 



E 



V V47T At 



exp {-[A[ -A- AtF q/ (6)] 2 /4At} 



In the Monte Carlo simulation, each iteration corresponds to a time step At. At each 
iteration, we sweep through each link in turn, and simulate the corresponding exponential 
factor in the sum on the right of ( p"2| ) by a random displacement of the link variable for each 
walker: 

AAi — At Fqi (8) + \, (12) 



where x is randomly chosen from a Gaussian distribution with standard deviation (2 At) . 
The first term in ( p"2| ) is the "drift" term, and the second is the "diffusion" term. The first 
exponential on the right of ( |l2| ) is simulated by simply multiplying the "weight" of each 
walker Wi by an equivalent amount. 

At the end of each iteration, the trial energy E^ is adjusted to compensate for any change 
in the total weight of all walkers in the ensemble; and a "branching" process is carried out, so 
that walkers with weight greater than (say) 2 are split into two new walkers, while any two 
walkers with weight less than (say) 1 /2 are combined into one, chosen randomly according to 
weight from the originals. This procedure of "Runge smoothing" ]22] | maximizes statistical 
accuracy by keeping the weights of all the walkers within fixed bounds, while mimimizing 
any fluctuations in the total weight due to the branching process. 

When equilibrium is reached after many sweeps through the lattice, the average value 
of the trial energy Et will give an estimate of the ground-state energy Eq, and the weight 
density of the ensemble in configuration space will be proportional to C'o^E't- Various cor- 
rections due to the finite time interval At have been ignored in this discussion, and the 
limit At — * oo must be taken in some fashion to eliminate such corrections. 

In the simulations presented here, a trial function for the ground state was chosen as 



-0t(6) = exp 



c 2^ cos Op + d 2_, cos(8p + 8p>) 
. p <pp'> 



(13) 



with two variational parameters c and d, where the sum over (PP 1 ) denotes a sum over 
nearest-neighbour pairs of plaquettes, forming rectangles (1x2 Wilson loops) on the lattice. 
That is, the uncorrelated single plaquette trial function used in |]l3|] has been extended to 
include a correlated term. In the single plaquette case (d — 0) it is straightforward to 
determine the optimal value for c [ p. 5| . For the correlated case (d ^ 0), the optimum 
values for c and d must be found by means of a separate Variational Monte Carlo (VMC) 
calculation, where the variational energy 

g vMcM)= /0M e) )a ge (14) 

is minimised. The procedure is described in |lq ] where a six- parameter variational wave 
function was used to estimate the mass gaps in the U(l) model at the VMC level. In 
short, one samples the local trial energy over configurations {0} randomly generated by a 
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Metropolis procedure, distributed with respect to the weight function (\I>t(@)) 2 - A down- 
hill simplex method was used to perform the minimisation. To this end -Evmc(c, d) was 
evaluated in a region centered around a guessed minimum (co,do) using the reweighting 
procedure described in fl5| , whereby uncorrelated configurations generated with respect to 
the weight function (\E f (©) c=Co ,^ =c / ) 2 were reweighted to give a distribution corresponding 
to the required weight function for a pair (c, d) in the neighbourhood of (co, do). This pro- 
cedure was iterated until the minimum (co,do) converged. We used 100,000 independent 
configurations on an 8 x 8 lattice to determine £Vmc(c, d). This allowed the optimal c and 
d to be fixed to within around 1%. Though it is possible to optimise c and d for each lattice 
size, for a given value of the coupling x, we used the same values of c and d for all lattice 
sizes L. The values of c and d used are listed in Table Q. 

The use of a correlated trial wave function is expected to markedly decrease the variance 
of estimators and the amount of branching in the population smoothing process, particu- 
larly in the physically interesting weak coupling regime (x large). Moreover, because the 
trial funct ion is closer to the true ground state, fewer forward walking iterations (see sec- 
tion I1IC) are required in order to derive converged estimators for correlation functions. 
These advantages come at a price in that updating of GFMC configurations becomes more 
complicated: The local trial energy is 

E L (8) = (4c- 2a;) ^ co s Op + 6<i ^ cos ( 6p + 9p ') ~ 4 F Ql> ( 15 ) 
p {PP') 1 

and the quantum force term is given by 

4 

Fqi = c (sin 9 Pl{ i } - sin6» Ql(i) ) + d^2 { sin ( 9 Pi(i) + e P 3 (i)) 

)}• (16) 



L (0 + P Q 3 W 



Here, P\(l) v ■ ■ ,Pa(1) and Qi(l),. ■ ■ ,Q&{1) denote the 4 closest plaquettes on either side of 
the link I, as illustrated in Fig. [l](b). These expressions are more difficult to compute than 
in the single parameter (d = 0) case. However, the increased complexity is easily offset by 
the gains from variance reduction. 



C. Forward- walking method 

The quickest method of estimating an expectation value (Q)o is simply to form the 
weighted average of Q over the ensemble of random walkers (we assume Q is diagonal 
in the chosen basis, for simplicity). This produces an estimate according to the distribu- 
tion 'J'r'fcb rather than $q. The estimate can be perturbatively improved ||, but there 
remains an unknown systematic error due to the dependence on the trial function v&t- F° r 
this reason, we have preferred to use the so-called "forward-walking" method for estimating 
expectation values. 

The forward-walking method is a robust technique for estimating expectation values 
pl| , p8|j39| , based on the following equation |2^] for an operator Q: 

<Q>a . (yay , 17 , 



(V T \K J Q\<i> ) 
)* T |*T J |$o) 

E M&j, ej-i) . . . k(e 2 , e 1 )g(e 1 )l»o(8 1 

E K(Qj, Qj-i) . . . K(Q 2 , 9i)$ (9i) 



(18) 
(19) 



where K (8 j, O is the evolution operator for time At, and K(Qj,Qj-i) is the same 
operator in the similarity transformed basis. Again we have assumed that the operator Q 
is diagonal in the basis of plaquette variables 8. 
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This equation is implemented by j2l| . |38| , |39j| the following procedure: 

i) Starting from the trial state, iterate until equilibrium is achieved, then begin a mea- 
surement; 

ii) Record the value <3(0;) for each walker ("ancestor") at the beginning of the measure- 
ment; 

iii) Propagate the ensemble as normal for J iterations, keeping a record of the "ancestor" 
of each walker in the current population; 

iv) Take the weighted average of the Q(0i) with respect to the weights of the descendants 
of Qi after the J iterations, using sufficient iterations J that the estimate reaches a 
'plateau'. 

This procedure has been tested for the magnetization in the 2D lattice Heisenberg model 
by RungeJ22j , and for correlation functions in the ID transverse Ising model by Samaras and 
Hamer |2Cj] , and works very well. The drawback to the procedure is that after a large number 
of iterations J many of the "ancestors" will die out, leaving no descendants, which leads to 
a progressive loss of statistical accuracy. Thus it is even more crucial in this connection to 
use a good guiding wavefunction 'J't so as to minimize "branching" . 



D. Timelike Correlations 



In order to estimate mass gaps in this model, we have again used the forward-walking 
technique to measure correlations between operators at different times. The mass gaps can 
then be estimated from the decay constants for these correlation functions. 

A timelike correlator (Qi(t)Q2(0))o can be found from the following equations: 

< Qi{t)Q2(v) >o= ^ [(f> (20) 



~ < ^ T \K J Q 1 K N Q 2 \^o > 

<* T |^ J+7V |$o> [ ' 

= ZK J Q 1 K N Q 2 <i>o m) 

where NAt = r. This can be implemented in much the same way as an expectation 
value (assuming both Qi and Q 2 are diagonal in the weak-coupling representation) . At the 



beginning of the measurement, record the 'ancestor' configurations as in Sec. (IIIC). Then 
allow the ensemble to propagate for time r, with the branching process turned off so that 
each state retains its identity. At the end of time r, record the initial value Q\ and the 
final value Q 2 . Propagate each state for a further J iterations as before, and then average, 
weighting each 'ancestor' state according to the forward- walking prescription of Sec. ( III C| ) . 



IV. RESULTS 



Simulations were carried out for L x L lattices up to L = 16 sites, using runs of 10,000 
walkers over 50,000 iterations. At each iteration several sweeps of the lattice were performed, 
after which Runge smoothing (the branching (combining) of high (low) weight walkers) was 
imposed on the walker population. The average percentage of walkers branched/combined 
per iteration depends on a number of factors: clearly, the time step At and number of 
lattice sweeps performed per iteration can be adjusted in order to control the extent of 
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the branching. We wish to choose values of Ar which are sufficiently small that time dis- 
cretization errors can be made negligible. Furthermore, the number of sweeps performed per 
iteration was chosen so that essentially uncorrelated measurements could be made roughly 
every 150 iterations. We found that, for the lattice sizes L and couplings x considered, 
this could be achieved in such a way that on average no more than 10% of the walkers 
were branched/combined per iteration. However, the extent of the smoothing per iteration 
depends on the quality of the variational guiding function. Generally, the two-parameter 
guiding function ( |l"3| ) works better for smaller couplings x and smaller lattice sizes L. Under 
the most strained conditions considered (x = 4, L = 16) around 10% of the walkers were 
branched/combined per iteration for the values of Ar chosen and the number of sweeps 
per iteration required to get independent measurements, whereas for x < 1 the ratio was 
a fraction of a percent. The first 2000 iterations were discarded to allow for equilibration. 
Two values of the time step Ar were used, differing by a factor of 5 (e.g. Ar = 0.02 and 
Ar = 0.005). The number of sweeps performed for the smaller Ar value was 5 times as 
large as that for the larger value and as a result measurements were equally uncorrelated 
in the two cases. The typical number of sweeps used per iteration, e.g. for Ar = 0.02 and 
Ar = 0.005 was 1 and 5 respectively. Linear extrapolations were made to the limit Ar = 0. 
The linear dependence on Ar has been checked previously jl3| . 



A. Ground-state Energy 



Results for the ground-state energy per site Eg / L 2 at coupling x = 4 are graphed against 
1/L 3 in Fig. ||. They agree within errors with those of the earlier study |l3| up to L = 10, 
and extend them to L — 16. It can be seen that the data are fitted quite well by the form 



E /L 2 



-4.414(1) 



2.53 



(23) 



which compares very well with effective Lagrangian theory and the weak-coupling series 
prediction pMl 



E /L 2 = -4.413(2) 



2.48(1) 



(24) 



The estimated bulk value —4.414(1) compares well with estimates of —4.43(2) from strong- 
coupling series |p8|, —4.415(6) from a ^-expansion p9(], and —4.412 from the coupled cluster 
method (CCM) 



B. Wilson Loops 

The expectation values of Wilson loops were computed using the forward-walking method 



outlined in section 1IIC. Measurements of the observables were made in cycles starting 
around every 150 iterations, with typically 10 forward-walking weighted averages over the 
ensemble taken at time steps of 12-15 iterations. For each time step the weighted averages 
were then block averaged over successive measurement cycles. Fig. ||(a) shows the forward- 
walking convergence of the 3x3 Wilson loop W(3, 3) measured on the 16 x 16 lattice 
at coupling x = 4 for a trial run involving 1200 walkers, 20,000 iterations, a time step 
Ar = 0.005, 3 sweeps of the lattice per iteration, with 12 forward-walking weighted averages 
being taken every 20 iterations in measurement cycles made 250 iterations apart. Results 
are shown for two different guiding functions, using the 1-parameter and 2-parameter forms 
respectively. It can be seen that in both cases the data relax exponentially towards a 
common equilibrium value, which can be estimated by making an exponential fit to the 
data. The equilibrium value is then taken as the final result for the Wilson loop. It can 
also be seen that the 2-parameter form reduces the variance and produces much more rapid 
convergence to the asymptotic value. 
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Typical results from a production run (using the 2-parameter guiding function) are shown 
in Fig. ||(b). In this case L = 12, x = 2 and 50,000 iterations are performed with f 0,000 
walkers. Measurement cycles involving f forward-walking weighted averages 16 iterations 
apart were started every 165 iterations. Results for two values of At (0.05 and 0.01), using 
1 and 5 lattice sweeps per iterations respectively, are shown. 

Fig. [| displays the values for the 'mean plaquette' P — W(l, 1) at coupling x = 4, graphed 
against 1/L 3 . There is evidently a discrepancy here between our present results and those 
of obtained using the secondary amplitude technique. Either the errors in Q have 
been underestimated for the larger lattices and couplings, or the discrepancy could have 
arisen from bias in the secondary amplitude technique Jl9| , p0[ . Our present results show a 
consistent finite-size scaling behaviour for P{L): 

22 

P{L) = 0.7584(4) + ^, (25) 
to be compared with the weak-coupling series prediction p3 

P(L) = -~[E /L 2 ] (26) 

= 0.7593(3) + ^^, (27) 

obtained using the Hcllmann-Fcynman theorem. The agreement is once again quite good. 
Extrapolation to the axis gives a bulk limit of P = 0.7584(4), to be compared with strong- 
coupling series §§ 0.80(3), the i-expansion §§ 0.757(4), and the CCM § 0.7585. The 
results for other Wilson loops scale similarly with lattice size. For example, the finite-size 
scaling of the 4x4 Wilson loop is shown in Fig. ||. 

Of course, because the model has a gap for any finite value of x, one should expect to 
see a crossover from this algebraic (free-photon theory) scaling to exponential convergence 
for large enough L. Indeed, for smaller values of the coupling, exponential scaling sets in 
rapidly enough that essentially bulk results can be obtained on relatively small lattices. This 
is illustrated in Fig. |^ in the case x = 2 where we see a definite crossover from algebraic 
to exponential scaling for the mean plaquette and the 2x2 Wilson loop when L reaches 
around 10. 

Table || lists the final estimates we have obtained for the bulk ground-state energy per 
site and Wilson loop values at some selected couplings. 



C. String Tension 



In a confining model, the Wilson loops are expected to behave asymptotically as 

W{m, n) ~ exp [— KA] 



(28) 



where A = mn is the area of the loop and K is the string tension. In Fig. [7] we illustrate 
this behaviour for the x = 2 case on the 12 x 12 lattice by plotting the Wilson loops W(n, n) 
and W(n,n — 1) as a function of area. The standard estimator of the string tension is the 
Creutz ratio 



K ~ R„ = - In 



W{n,n)W(n- l,n- 1) 



W(n,n- l) 2 



(29) 



Fig. ||(a) shows the Creutz ratios R n graphed against 1/ A (where A = n(n— 1) is the average 
area of the loops used to form the ratio) at x = 4, from which we obtain the rough estimate 



K = 0.05(2) 



(30) 



Because the Creutz ratios decrease very rapidly with n, and the relative error correspond- 
ingly increases, it is very difficult to determine a reliable value for the string tension. Given 
the limited number of data points, it is probably more fruitful to use a two-point estimator 
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K-K = --ln\ ^ n \ ] (31) 
n \W{n,n - 1)J 

formed from successive pairs of Wilson loops. Included in Fig. |§|(a) are plots of these 
estimators against 1/ A. Unfortunately, the two-point estimators display oscillations (which 
the Creutz ratios were designed to eliminate) so it is again difficult to perform a linear 
extrapolation to obtain an estimate of the string tension. Nevertheless, a linear fit through 
the points gives an estimate of the string tension of K sa 0.03, somewhat lower than the 
result above. 

We have performed a similar analysis for x = 2, working with the 12 x 12 lattice which is 
sufficiently close to the bulk limit. The results for the string tension estimators are shown in 
Fig. |^(b). In this case a linear extrapolation of the two-point estimators yields an estimate 
of K rj 0.08 for the string tension. 

The string tension has previously been calculated as an energy per unit length, a. The best 
available values at x = 2 are a = 0.28(1) from stochastic truncation jyjj, and a = 0.282(2) 
from an 'exact linked cluster expansion' p0[ . The quantities are related by 

a = vK (32) 

where v is the speed of light, estimated as v — 2.27 at a; = 2 from the weak-coupling 
expansion ^4j. Hence a — 0.28(1) corresponds to K = 0.12(1). This is in rough agreement 
with the result above. 

We also attempted the analysis for smaller values of x. However, because the Wilson 
loops decay extremely rapidly in these cases, it is not possible to obtain enough data points 
to make sensible extrapolations. 



D. Mass Gaps 

We have attempted to estimate mass gaps from the exponential decay of correlation 
functions, as is done in the Euclidean approach. The lowest-lying excitations in the strong- 
coupling limit are the single-plaquette excitations, antisymmetric or symmetric under reflec- 
tions respectively. At first we attempted to use the spatial correlations between plaquette 
operators for this purpose; but the signal from the connected spatial correlations turns out 
to be very small, and can even change sign, giving no good exponential decay signal. We can 
understand this problem as follows: : we are restricted to a plane in the spatial directions, 
and can therefore only measure correlations between plaquette operators "edge on" to each 
other, as illustrated in Fig. |^(a), as opposed to the "face on" configuration in Fig. ^(b). 
This means there is very poor overlap with the configurations of interest in the intermediate 
state. 

We are thus compelled to fall back on the correlations in imaginary ti me t, which do 



correspond to "face on" plaquettes as in Fig. |9j(b). Using the method of |IID|, we have 
measured correlation functions 

/w = <(E p (°))(E p '( r )))o ( 33 ) 
p pi 

where the plaquette operators have been summed over all positions to project out the zero- 
momentum intermediate states, and the sums are either symmetric or antisymmetric under 
reflections. 

As with the Wilson loop calculations, typically 50,000 iterations were performed with 
10,000 walkers and forward-walking measurment cycles were started every 150 iterations. 
As mentioned, at the start of each measurement cycle, Runge smoothing was switched 
off so that no walkers were created or destroyed while timelike correlation functions were 
calculated. Between 10 and 20 timelike correlations (between the plaquette operators at the 
start of the measurement cycle r = and later times r)were measured at intervals of one 
or two iterations. For each of the timelike measurements the "physical time" r is given by: 

r = At x Asweep x Ainterval 
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where At is the basic timestep, iVswEEP is the number of sweeps of the lattice and 
•^interval denotes the number of iterations between the calculation of initial and final 
plaquette operators. 

At this point Runge smoothing was switched on again and forward-walking weighted 
averages of the timelike correlations were calculated over most of the measurement cycle 
at intervals of 10-20 iterations. Again the weighted averages were block averaged over all 
measurement cycles, and exponential fitting was used to find the forward-walking relaxed 
estimates of the timelike correlation functions. As with the Wilson loops, two values of the 
time step At were used and linear extrapolation was used to obtain At = estimates of 
the correlators. Again, the two values of At used differed by a factor of 5, and 5 times as 
many lattice sweeps were performed per iteration for the smaller At so that the "physical" 
times t over which the correlations /(t) were measured were the same (i.e. At x A*sweep 
was kept fixed). 

A typical plot of the resulting estimates for the timelike correlation function /(t), as a 
function of r, is shown in Fig. [H] for the antisymmetric correlator with x = 0.5 and L = 4. 
It can be seen that /(t) shows a smooth exponential decay with t. Defining an 'effective 
mass' for each pair of points by 

m(r) =-1u[/(t)//(t-1)], (34) 

one finds that the effective mass decreases somewhat with r, as might be expected. An ad 



hoc fit in 1/t, illustrated in Fig. 11, gives a final estimate for the mass of the corresponding 



intermediate state. These values are listed as a function of lattice size in Table II] for 
couplings x = 0.5, 1 and 2. For the higher couplings the exponential decay is so slow 
that we have been unable to obtain useful results for the mass. To do so would require 
measurements to be made over substantially larger time scales with significantly greater 
accuracy. Unfortunately, our results were not sufficiently accurate to sensibly resolve the 



form of the lattice size dependence for the couplings studied. In Table |TJ we list results 
of estimates of the bulk mass from other sources. It can be seen that for the couplings 
considered, the agreement between methods is reasonable. 



V. SUMMARY AND CONCLUSIONS 



In this paper we have shown how the standard 'forward-walking' methods of quantum 
many-body theory j^] can be used to calculate expectation values and correlation functions in 
Hamiltonian lattice gauge theory. Accurate values for the Wilson loops have been obtained, 
and their finite-size scaling behaviour has been demonstrated. The string tension can then 
be estimated using two-point and Creutz ratios. The Wilson loop values dropped away 
too quickly with lattice size to give accurate string tensions at small couplings x; but this 
problem could presumably be rectified by looking at timelike correlators — either timelike 
Wilson loops, or timelike correlations between Polyakov loops. 

We have also shown how mass gaps can be estimated from the exponential decay of time- 
like correlators. This approach is very reminiscent of the techniques used in the Euclidean 
regime: our freedom of choice of the timelike measurement interval is akin to the choice of 
an anisotropic lattice in the 'timelike' direction in the Euclidean case. Reasonable results 
have been obtained for the masses at small couplings x, but not at large x. 

The accuracy of the final results for the string tension and mass gaps is not remarkable; 
but it must be recalled that this is the first time that unbiased estimates have been obtained 
for these quantities using a 'weak-coupling' representation and GFMC techniques. There 
are no doubt many ways by which the results could be improved: for example, by the use of 
improved guiding wavefunctions, or improved ("fuzzed") correlation functions, or improved 
Hamiltonian operators. Our aim here has been to 'prove the concept', that forward-walking 
techniques can give reliable estimates for the correlation functions. In future work we plan 
to apply similar methods to the (3+l)D SU(3) Yang-Mills theory. 
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x c d 

Ol 0.239 0.0088 

1 0.409 0.033 

2 0.585 0.089 
4 0.771 0.150 



TABLE I. Optimal values ol the coefficients c and d used in the guiding lunction (Jl 3|) - 

X = 1 



L = 4 L = 6 L = 8 

£0/^ -0.46520(3) -0.465127(2) -0.465065(2) 

W(l,l) 0.4337(1) 0.43350(7) 0.43346(6) 

W(l,2) 0.2161(1) 0.21526(7) 0.21521(6) 

W(2,2) 0.0649(1) 0.06377(7) 0.06369(6) 

W(2,3) 0.01946(5) 0.01927(4) 

W(3,3) 0.00391(5) 0.00378(5) 

W(3,4) 0.00077(3) 



x = 2 



L = 4 L = 6 L = 8 £ = 10 L = 12 

£ /£ 2 -1.58249(9) -1.5707(6) -1.56993(5) -1.56996(3) -1.56998(3) 

W(l,l) 0.6492(3) 0.6394(2) 0.6378(2) 0.6380(2) 0.6378(2) 

W(l,2) 0.4772(3) 0.4572(2) 0.4544(2) 0.4545(2) 0.4545(2) 

W(2,2) 0.3134(3) 0.2726(2) 0.2675(2) 0.2673(2) 0.2675(2) 

W(2,3) 0.1683(2) 0.1623(2) 0.1616(2) 0.1615(2) 

W(3,3) 0.0932(3) 0.0855(2) 0.0846(2) 0.0837(2) 

W(3,4) 0.0459(2) 0.0449(2) 0.0440(2) 

W(4,4) 0.0218(2) 0.0208(2) 0.0205(2) 

W(4,5) 0.0097(2) 0.0099(3) 

W(5,5) 0.0036(4) 0.0041(5) 



x = 4 



L = 4 L = 6 L = 8 L = 10 L = 12 L = 16 L = oo 

Bo/L^ -4.4545(1) -4.4258(2) -4.4189(2) -4.41653(8) -4.4157(1) -4.4146(6) -4.414(1) 

W(l,l) 0.7624(2) 0.7592(3) 0.7588(2) 0.7587(2) 0.7586(4) 0.7588(6) 0.7584(1) 

W(l,2) 0.6338(4) 0.6256(4) 0.6236(2) 0.6234(4) 0.6223(4) 0.6233(7) 0.6226(2) 

W(2,2) 0.4980(5) 0.4751(7) 0.4706(4) 0.4693(6) 0.4660(6) 0.4659(9) 0.4665(2) 

W(2,3) ■ 0.3739(8) 0.3651(4) 0.3622(7) 0.3575(6) 0.356(1) 0.3573(3) 

W(3,3) 0.2872(9) 0.2713(5) 0.2662(8) 0.2598(8) 0.257(1) 0.2581(5) 

W(3,4) 0.2066(5) 0.1981(8) 0.1925(9) 0.189(1) 0.1883(7) 

W(4,4) 0.1552(6) 0.1426(9) 0.138(1) 0.133(2) 0.130(1) 

W{4,5) 0.1051(9) 0.101(1) 0.094(2) 0.091(1) 

W(5,5) 0.076(1) 0.071(1) 0.066(2) 0.064(2) 

W(5,6) 0.051(1) 0.045(2) 0.040(3) 

W(6,6) 0.034(1) 0.031(3) 0.029(5) 



TABLE II. Estimates of the ground-state energy per site Eo / L 2 and Wilson loop values W(m, n) 
as functions of lattice size L for various couplings x. In the x = 4 case we extrapolate to the bulk 
limit L — 00 by fitting to the form Wz,(m, n) ~ ^^(m,?!) + A/L 3 . 
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FIG. 1. (a) The links h,. . . ,U associated with a plaquette P on the lattice, (b) The plaquettes 
P\{V),. . . ,Pa(1) and Qi(l),. . . ,Qa(1) associated with a link I on the lattice. 
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FIG. 2. Ground-state energy per site graphed against 1/L 3 , where L is the lattice size, at 
x — 4. The triangles show our present estimates (GFMC), the diamonds are from the stochastic 
truncation calculations from ref. The solid and dashed lines are least squares fits to the 

GFMC and stochastic truncation results respectively. The dotted line is from the predictions of 
weak coupling theory (see text). 
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FIG. 3. (a) Forward- walking weighted averages of the Wilson loop W(3, 3) at coupling x = 4 and 
lattice size L = 16, as a function of the number of forward-walking steps. Each forward walking step 
involves 20 MC iterations which in turn consist of 3 lattice sweeps, with a time step of At — 0.005. 
Open diamonds are results obtained from a 1-parameter guiding function, filled diamonds from a 
2-parameter guiding function, (b) Convergence of forward-walking weighted averages for the 3x3 
Wilson loop for x = 2 on the 12 x 12 lattice. Each forward walking step consists of 16 MC iterations. 
The open diamonds are for At — 0.05 (1 lattice sweep per MC iteration) and the closed diamonds 
are for At = 0.01 (5 lattice sweeps per MC iteration). 
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FIG. 4. The mean plaquette value W(l, 1) as a function of 1/L 3 , at x — 4. The triangles are 
our present forward-walking GFMC results; the diamonds are GFMC results from ref. [131 using 
the "secondary amplitude" method of evaluating expectation values. The solid and dashed lines are 
least squares fits to the forward-walking and secondary amplitude results respectively. The dotted 
line is from the predictions of weak coupling theory (see text). 




1G 



0.155 
0.15 
^ 0.145 
0.14 

0.135 



X = 4 



A' 



0.0005 



0.001 

1/L 3 



0.0015 



0.002 



FIG. 5. Estimates of the 4x4 Wilson loop W{A, 4) at x = 4, as a function of 1/L 3 . The dashed 
line is a least squares linear fit to the data. 
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FIG. 6. (a) The mean plaquette, W(l, 1) and (b) the 2x2 Wilson loop as functions of 1/L 3 for 
the x = 2 case. 
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FIG. 7. Wilson loops W(n, n) and W(n,n — 1) as functions of their area for x = 2 on the 12 x 12 
lattice. 



0.1 
0.09 
§ 0.08 
§ 0.07 
a 0.06 

c 
o 

S 0.05 

a 

u 

" 0.04 



so 
c 

'fi 



0.03 



0.02 



(a) x = 4 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

1/A 



0.2 
0.18 
2 0.16 
.1 0.14 



c 



1) 



0.12 



0.1 



bO 08 
c 

'£ 0.06 

GO 

0.04 E 



(b) x = 2 



0.1 



0.2 



0.3 



0.4 



0.5 



0.8 



0.6 0.7 
1/A 

FIG. 8. 'Two-point' string tension estimates R n i (diamonds) and Creutz ratios R n (triangles) 
for (a) x = 4 (using extrapolated L = oo estimates of the Wilson loops) and (b) x = 2 (using 
L = 12 estimates of the Wilson loops), graphed against 1/A, where A is the average area of the 
Wilson loops used to form the ratios. 
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FIG. 9. (a) Two plaquettes "edge on" on a lattice with two spatial dimensions, (b) "Face on" 
plaquettes on a lattice with three spatial dimensions. 
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FIG. 10. Timelike correlation function /(r) (for antisymmetric states) as a function of imaginary 
time r, for coupling x — 0.5 and lattice size L = 4. 
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FIG. 11. "Effective mass" ms(r) for the symmetric state as a function of 1/r for coupling 
x = 0.5 and lattice size L — 4. (Here r is measured in terms of MC iterations which correspond to 
a physical time of 0.15). 
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